home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Languages / RLaB 1.18c / testmatrix / dorr.r < prev    next >
Text File  |  1994-12-20  |  2KB  |  61 lines
  1. //-------------------------------------------------------------------//
  2.  
  3. // Synopsis:    Dorr matrix - diagonally dominant, ill-conditioned, 
  4. //        tridiagonal.
  5.  
  6. // Syntax:    DL = dorr ( N , THETA )
  7.  
  8. // Description:
  9.  
  10. //    dorr(N, THETA) returns a list containing the vectors defining
  11. //    a row diagonally dominant, tridiagonal M-matrix that is
  12. //    ill-conditioned for small values of the parameter THETA >= 0.
  13. //      The columns of INV(C) vary greatly in norm.  THETA defaults to
  14. //    0.01. The amount of diagonal dominance is given by (ignoring
  15. //    rounding errors): 
  16.  
  17. //        COMP(C)*ONES(N,1) = THETA*(N+1)^2 * [1,0,0,... 0,1]'.
  18.  
  19. //      Reference:
  20. //      F.W. Dorr, An example of ill-conditioning in the numerical
  21. //      solution of singular perturbation problems, Math. Comp., 25 (1971),
  22. //      pp. 271-283.
  23.  
  24. //    This file is a translation of dorr.m from version 2.0 of
  25. //    "The Test Matrix Toolbox for Matlab", described in Numerical
  26. //    Analysis Report No. 237, December 1993, by N. J. Higham.
  27.  
  28. //-------------------------------------------------------------------//
  29.  
  30. dorr = function ( n , theta )
  31. {
  32.   local (n, theta)
  33.  
  34.   if (!exist (theta)) { theta = 0.01; }
  35.  
  36.   c = zeros(n,1);
  37.   e = c;
  38.   d = c;
  39.  
  40.   // All length n for convenience.  Make c, e of length n-1 later.
  41.  
  42.   h = 1/(n+1);
  43.   m = floor( (n+1)/2 );
  44.   term = theta/h^2;
  45.  
  46.   i = (1:m)';
  47.     c[i] = -term*ones(m,1);
  48.     e[i] = c[i] - (0.5-i*h)/h;
  49.     d[i] = -(c[i] + e[i]);
  50.  
  51.   i = (m+1:n)';
  52.     e[i] = -term*ones(n-m,1);
  53.     c[i] = e[i] + (0.5-i*h)/h;
  54.     d[i] = -(c[i] + e[i]);
  55.  
  56.   c = c[2:n];
  57.   e = e[1:n-1];
  58.  
  59.   return << c = c; d = d; e = e >>;
  60. };
  61.